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ABSTRACT 

The scale-invariant glitch statistics observed in individual pulsars (exponential 
waiting-time and power-law size distributions) are consistent with a critical self- 
organization process, wherein superfluid vortices pin metastably in macroscopic 
domains and unpin collectively via nearest-neighbor avalanches. Macroscopic in- 
homogeneity emerges naturally if pinning occurs at crustal faults. If, instead, 
pinning occurs at lattice sites and defects, which are macroscopically homoge- 
neous, we show that an alternative, noncritical self-organization process operates, 
termed coherent noise, wherein the global Magnus force acts uniformly on vor- 
tices trapped in a range of pinning potentials and undergoing thermal creep. 
It is found that vortices again unpin collectively, but not via nearest-neighbor 
avalanches, and that, counterintuitively, the resulting glitch sizes are scale in- 
variant, in accord with observational data. A mean- field analytic theory of the 
coherent noise process, supported by Monte-Carlo simulations, yields a power- 
law size distribution, between the smallest and largest glitch, with exponent a in 
the range —2 < a < 0. When the theory is fitted to data from the nine most ac- 
tive pulsars, including the two quasiperiodic glitchers PSR J0537— 6910 and PSR 
J0835— 4510, it directly constrains the distribution of pinning potentials in the 
star, leading to two conclusions: (i) the potentials are broadly distributed, with 
the mean comparable to the standard deviation; and (ii) the mean potential de- 
creases with characteristic age. Fitting the theory to the data also constrains the 
pinned vortex fraction and the rate of thermal creep. An observational test is pro- 
posed to discriminate between nearest-neighbor avalanches and coherent noise: 
the latter process predicts a statistical excess of large glitches ('aftershocks') fol- 
lowing a large glitch, whereas the former process does not. Its discriminatory 
power is discussed under various microphysical scenarios. 

Subject headings: dense matter — hydrodynamics — stars: interior — stars: 
neutron — stars: rotation 
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Introduction 



The number of pulsar glitches recorded historically has quadrupled in the last five years 
and now approaches 300 events, following improvements in the sensitivity and duty cycle of 
radio pulsar timing, pioneered b y the Parkes Multibeam Su r vey jHobbslbooil iKrawczvk et al.l 



2006 



Per alt a 



2003 : [Manchester et al.l l2005l : iJanssen fc Stapperj l2006l : iMiddleditch et al 
20071 ). It is now possible to disaggregate the data and measure reliably the glitch size 



and waiting-time distributions in individual pulsars. iMelatos et al.l (120081 ) showed that the 
distributions are consistent with power laws and exponentials respectively in nine objects 
(two of which also feature a smaller quasiperiodic component), suggesti ng that the glitch 



mechanism is predorn i nantly scale invarian t and obeys Poisson statistics (jWong et al. 
Alpar fc Bavka3l200d : iMelatos et aPboOsh . 



2001 



Scale invariance and Poisson statistics are universal characteristics of self-organized 
critical systems, in which discrete, interacting elements adjust in response to a slow, local- 
ized, external driver via intermittent avalanches, i.e. n earest-neighbor "domino cha ins" of 
local, impulsive, threshold-activated relaxation events (jJenseru Il998l : ISornettd l2004f) . Self - 
organized cri ticality is obser ved widely in nature, for example in sandpiles fBak et al.lll987n . 
earthquakes ( Sornette 2004), sola r fiares (Lu fc Hamilton 1991 : Wheatland boOO ) . and mag- 
netized type II superco nductors (IField et al.l 119951). Two popular pulsar glitch paradigms. 



involving crust fracture (lAlpar et al. 



of superfiuid vortices (lAnderson &: Itoh 



1996 



naturally to avala nche dynamics (IMe. 



cellular automata (IMorley &: Schmidt 



Middleditch et al. 



19751: 



atos et al. 



1996 



Cheng et al 



2006h and collective unp inning 



19881 : lAlpar et al.lll996h . lead 



20081) and can be simula ted efficiently with 



Warszawski fc Melatod 120081 ) . 



As avalanches traverse a self-organized critical system, they leave behind inhomogeneities 
on all scales, up to and including the system size. In the vortex unpinning paradigm for pulsar 



(Cheng et al. 


1988: 


Alpar et al. 


1996: 


Wong et al. 


2001) 



voirs occupy a sizable portion of the star. This picture works admirably if pinning occurs 
along a macroscopic netw ork of crustal faults (created by seismic disturbances, for example) 



(IMiddleditch et al.ll2006l ). It is harder to reconcile with pinning at nuclear lattice defects 



(e.g. interstitial vacancies) or simply lattice sites, which are separated by tens of lattice spac- 



ings and are therefore homo geneous on macroscopic scales (|Joneslll998l : lDonati fc Pizzochero 



20031 : lAvogadro et al.l 120071 ). In this picture vortices hop between adjacent pinning sites in 
response to thermal fiuctuations and the Magnus stress (from differential spin down of the 
crust and superfi uid) without deviating macroscopically from a regular Abrikosov lattice 



(iLink et al.lll993f ). 



A central puzzle concerning vortex unpinning is how it accounts for the scale invari- 
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ance of pulsar glitches. Normally, scale invariance is a sign of collective behavior, involving 
correlations between inhomogeneities on all scales, including the largest. But, as explained 
above, large-scale inhomogeneity is hard to contrive when pinning occurs at lattice defects or 
lattice sites, so the inter-vortex forces approximately balance. How, then, do vortices unpin 
in sympathy in such large (macroscopic) numbers? Furthermore, why do glitch sizes vary by 
many orders of magnitude (four decades in one object) from glitch to glitch? We emphasize 
at this point that the size of a glitch is related to not only the number o f vortices that un- 



pin but also the distance moved in the radial direction before repinning (lAlpar et al.lll984 



Cheng et al.lll988l : I Jahan-Miril |2006| ) . If nuclear lattice pinning is characterized by a t ypical 
pinni i ig potentiator even a moderate range of pinning potentials, as theory predicts (jJones 
19981 : iDonati fc Pizzocherd l2006l : lAvogadro et al.l 120071 ) , and if the global Magnus stress is 
felt uniformly by all pinned vortices, one naively expects glitches in any individual pulsar 
to recur periodically (whenever the Magnus stress rises to match the typical pinning force) 
and hence display approximately equal sizes, contrary to what is seen. These are profound 
challenges for all glitch mechanisms. 



In this paper, we invoke the coherent noise process introduced by ISneppen &: Newman 



( 119971 ) to describe the collective dynamics of vortex unpinning and calculate the resulting 
glitch statistics. Remarkably, we find that scale invariance emerges automatically between 
the minimum and maximum sizes, even without nearest-neighbor avalanches and large-scale 
inhomogeneity. In ^ we specify the model and calculate the glitch size distribution ana- 
lytically from first principles, in the stationary, mean-field approximation. The results are 
expressed in terms of three principal variables: the mean glitch rate, which is observable; the 
mean thermal unpinning rate, which does not equal the mean glitch rate; and the distribu- 
tion of pinning potentials, which can be predicted from nuclear physics. In ^ we compare 
the analytic theory against the output of Monte-Carlo simulations and show how existing 
and future observational data can be used to constrain the distribution of pinning potentials. 
By way of illustration, the theory is fitted to data from the nine most active glitchers cur- 
rently known. In ^ we explore the implications of our results for the nuclear microphysics 
of vortex pinning. We also propose an observational test that may discriminate between 
coherent noise and nearest-neighbor avalanches. 



We emphasize, at the outset, that the coherent noise mechanism does not operate under 



^The term 'coherent noise' is something of a misnomer when apphed to vortex unpinning. The noisy 
element, namely thermal creep, does not operate coherently throughout the star; rather, the vortices feel 
a globally coherent Magnus stress. Arguably, the adjective 'coherent stress' describes the model better. 
However, we elect to retain the original terminology in this paper in order to preserve consistency with the 
statistical mechanics literature. 
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conditions when pinning does not occur (jjonea 119971 . Il998l ) 



Coherent noise mechanism 



Sneppen fc Newman! (1l997l ) first postulated the coherent noise mechanism to describe 
discrete, far-from-equihbrium systems driven by a global stress, which acts on all elements 
in the system simultaneously. Elements respond by relaxing locally via threshold-activated 
events. Coherent noise offers a pathway to self-organization which is neither critical nor 
interaction-dominated; for instance, the stationary state is homogeneous on average, without 
long-range spatial correlations mediated by nearest-neighbor interactions. By contrast, in 
critical systems hke sandpiles, the external driver is localize d, and its influence propagate s 
across the system over time via scale-invariant avalanches (IBak et al.l 119871 : iJensenI Il998l ). 
Nevertheless, counterintuitively, coherent noise does produce intermittent, collective events 
with a scale-invariant size distribution, just like self-organized critical systems. 

In this section, we describe a simple cellular automaton that models the collective dy- 
namics of vortex unpinning in a neutron star as a coherent noise process ( §2.11 and §2.2p . 
We then solve analytically for the mean-field behavior of the model ( §2.31) and derive the 
distribution of glitch sizes as a function of the pinning parameters ( §2.41) . 



2.1. Pinning at lattice sites and defects 

Consider a neutron star containing N pinned superfluid vortices, each carrying circula- 
tion K, amounting to a fraction e of all the vortices in the star. In this paper, we take N to be 
constant, as we seek to model the short-term glitching of individual pulsars, monitored over 
decades, rather than the long-term glitching of the whole pulsar population, which evolves 
on the spin-down time-scale z//z>. Here, u is the spin frequency. 

The microscopic cause of pinning, whether it be the attractive force of the lattice nuclei 
themselves, defects in the lattice, or macroscopic faults in the stellar crust, is an important 
factor in determining the degree of homogeneity of pinning sites. Pinning at lattice nuclei, 
the default assumption in this paper, i s the most homogeneo us. Interstitial pinning at the 



midpoint between nuclei (if it exists; see lAvogadro et al.ll2007l ) is the next most homogeneous 



^ A word of caution: self-organized criticality and coherent noise are not the only collective mechanisms 
that lead to intermittency and scale invariance in far-from-equi librium systern s. Other examples include 
percolation, multiplicative noise, and highly optimized tolerance (|Sornettdl2004) . 
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level of pinning. Actual defects (e.g. nuclear imp urities like monovacancies or shear layers) 
are estimated to occur roughly every ~ 30 nuclei (Ide Blasio &: Lazzarilll998l ). but these too 
are essentially "perfectly homogeneous" on the vortex lattice scale that matters in this paper. 
Macroscopic faults in the solid crust , e.g. "tectonic" plates formed by large-scale cracking (if 



it exists; see 



Middleditch et al.ll2006l ) are probably inhomogeneous on the vortex lattice scale 



and lie outside the scope of any coherent noise model. Suppose each pi nned vortex (labelled 
i, wit h 1 < i < A^) occupies the site of a defect in the nuclear lattice (jde Blasio fc Lazzari 
1998 ). whose pinning threshold (expressed as a force per unit length) is denoted by Fp\ 



If the Magnus force Fu exceeds Fp \ then the vortex unpins and moves with the local 
superfluid flow. Howe ver, pinning sites are abundant microscopically, occurring once every 
~ 30 lattice spacings (|Joneslll998l : iLink &: Cutlerll2002l ). Consequently, the unpinned vortex 
moves at most ~ 10~^^m, much less than the mean vortex separation (/€/47rz/)^/^, before it 
immediately repins at a new defect, with a new pinning threshold Fp^ Fp*^ in general]. The 
microscopic abundance of pinning sites is the key difference between this scenario and one 
driven by crust fracture, to which avalanche models apply; at all times, the pinned vortices 
form something close to a regular Abrikosov array, which, like the underlying pinning sites, 
is homogeneous on large scales. § 

Let (f){Fp)dFp be the fraction of defects whose pinning threshold lies between Fp and 
Fp + dFp. In principle, it is possible to predict (t>{Fp) ah initio from nuclear structure 
calculations, although attempts to do so have yielded conflicting results due to the subtlety of 
the physics. One approach, based on the local density approximation, suggests that pinning 
in the inner crust is strong est a t intermedia te densities and is predorn inantly interstitial 



(IDonati fc Pizzocherd l2003l . |2004J ) or nuclear (IDonati fc Pizzocherd l2006l ) . It yields pinning 
energies Ep in the range 1 MeV Fp < A MeV, although Ep can be ~ 20 times greater 
if the normal component of the superfluid is absent ("pure phase"). A second approach, 
based on the mean-field, Hartree-Fock-Bogoliubov approximation, su ggests that pinriing is 



strongest at low (5 x 10^^ kg m ^) and high (2 x lO^^kgm ^) densities ( iLink fc Epsteinlll991 



Avogadro et al.l 120071 ). It yields 1 MeV ^ -Ep ^ 3 MeV for Fermi momenta in the range 0.5- 
1 fm~^. Yet another line of argument suggests that pinning is too weak to occ ur at all, unless 
the concentration of monovacancies is unexpectedly high (j Jones! 119971 . 119981 ). Summarizing 
the foregoing results, we adopt as a crude working hypothesis the top-hat distribution 

0(Fp) = (2A)-iiJ(Fp - Fo + A)H{-Fp + Fq + A) , (1) 

where Fq is the mean, A is the half-width, and H{. . . ) is the Heaviside step function. 



•^In this paper, we declare the vortices to be rectihnear, although recent work suggests that the meridional 
flow induced by crust-superfluid differential rotation excites the Gla berson-Donnelly instability and generates 
a vortex tangle |Peralta et allboosl bood iMelatos fc PeraltalboOTl ). 
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(Note that we assu me Fp > in this paper.) Typical values are Fn = 1 x 10^^ N and 
A = 6 X lO^^Nm"^ (jPonati &: Pizzocherdl2006l : lAvogadro et al.ll2007l ). taking the superfiuid 
coherence length to be ^ = 5 fm. 

Remarkably, any physically reasonable form of (p{Fp) (e.g. top hat, power law, Gaussian) 
gener ates a scale invariant distribution of event sizes over some interval (jSneppen fc Newman 
19971 ). We explain why in ^ This robustness is an attractive feature of the coherent noise 
mechanism as a model for pulsar glitches. Note that ([T]), as written, is independent of 
position within the star. We make this simplification here in anticipation of the mean-field 
analysis to be carried out in §2.3[ but of course it is not realistic; pinning theories favor 



certain ranges of superfiuid density (jPonati &: Pizzocherdl2003l . 12004 l2006l : lAvogadro et al. 



20071 ). We will generalize ([T]) by letting it vary with position in a forthcoming paper. 



2.2. Cellular automaton for vortex unpinning 



Let us now evolve the pinned vortices in discrete time steps At, with the aid of a simple 
cellular automaton. At each time step, the following four rules are applied to update the 
state of the automaton. They encode the microphysics of vortex unpinning in an idealized 
fashion, but the resulting collective behavior is insensitive to th e details of the microphysics , 
as much accumulated experience with cellular automata shows (I Jensenlll998l : ISornettell2004j ) . 



1. A value of the global stress (here, the Magnus force) Fm is chosen at random from a 
probability distribution function iIj{Fm)- As the Magnus force originates from crust- 
superfiuid differential spin down, ip{Fu) shares the same form as the observed glitch 
waiting-time distribution (see below). 

2. A small fraction / -C 1 of the pinned vortices unpin at random, e.g. due to thermal 
fluctuations. For simplicity, we do not bias this random process towards particular 
pinning sites. A more realistic model might preferentially unpin those vortices with 
Fp*-* just above Fm, e.g. with Fm < Fp*^ < Fm + ksT/^'^, where T is the temperature 
of the crust and is Boltzmann's constant. 

3. The stress Fm is applied simultaneously to all the remaining {1 — f)N pinned vortices, 
viz. F^^ = Fm. All vortices with Fp*^ < Fm unpin. For simplicity, we take Fm to 
be uniform, but it is straightforward to let it vary realistically with distance from the 
rotation axis. Generalizing the model in this way does not alter its collective behavior 
at all. 
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4. Each unpinned vortex repins almost immediately at a nearby defect (see §2.ip and is 
assigned a new threshold. 



In the simplest version of the automaton, the pinned vortices do not interact, except 
implicitly through the mutual repulsion which keeps them in a regular Abrikosov array. 
Consequently, in the stationary state, long-range spatial correlations do not emerge, and 
the occupied pinning sites are distributed homogeneously (see ^ transitory correlations 
can arise by accident, of course). Macroscopic homogeneity is preserved even if nearest- 
neighbor vortex interactions are allowed, provided that the interactions are weaker than 
the global stress most of the time. On the other hand, in the interaction-dominated regime, 
homogeneity breaks down, ari d the coherent noise process transitions to an avalanche process 
jWarszawski fc Melatoslboosl ). 



What is the form of iI){Fm)'^ There is strong empirical evidence that glitch waiting times 



(Wone et al. 


2001; 


Melatos et al. 


2008: 



of anomalous braking indices in terms of unresolved glitches ( iJohnston fc Gallowayl Il999 



Alpar fc Baykalll2006l ). Hence the automaton is asynchronous; At, the time since the last 
glitch, is different at every time step and exponentially distributed. A physically reasonable, 
monotonically decreasing form of ip{Fy[) is therefore the Poisson-like distribution 



^(Fm) = 0- exp(-FM/a) 



(2) 



where a is a characteristic value of the Magnus stress that accumulates prior to a glitch. 
Possible definitions of a in terms of the physical parameters of neutron stars are discussed 
in § 14. 1[ For the remainder of this section and §121 cr is treated as a stress scale factor when 
analyzing the behaviour of the model, which depends on a only through the dimensionless 
combinations F^/a and A/ a. Crucially, the collective dynamics are the same irrespective 
of the exact form of ipiFm)- ISneppen fc NewmanI (119971 ) showed that any distribution that 
falls off sufficiently rapidly at large Fm, such that dxif^^x) scales as ['?/'(Fm)]" (a real) to 
leading order, generates a power-law distribution of event sizes over some size interval. 

In the absence of nearest-neighbor interactions, one must have f ^ 0. Otherwise, 
the system stagnates ultimately, with (Fu) <^ Fp'^ for all i and vortices unpinning ever 
more rarely as time passes. The competition between coherent forcing and thermal creep 
is therefore essential for the system to develop scale invariance. In more elaborate models 
incorporating nearest-neighbor interactions, thermal creep is strictly unnecessary and one 
can set f = 0, but, as discussed in ffl this is not believed to be the situation in pulsars 
( Link et al. 1993 ; Link fc Epstein 19961 ). Note that / must be small to account for a power- 
law distribution of event sizes over several decades (ISneppen &: NewmanI 119971 ). 
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2.3. Mean-field master equation 

We now analyse the cellular automaton in the mean-field approximation, exploiting the 
property of homogeneity. Let g{Fp)dFp be the time- averaged fraction of pinned vortices 



trapped in pinning sites with thresholds in the range Fp < F^'^ < Fp + dFp. Crucially, g{F^ 



is not the same as (f){Fp); the latter function counts all pinning sites without bias, whereas 
g{Fp) counts only those sites that are actually occupied in the stationary state. To illustrate 
this key point, consider two sites i and j, with Fp^ = Fq — A/2 and Fp^ = Fq + A/2. If 
4>{Fp) is given by (Q, say, then we have ^f-Fp*"*] = 0[Fp"'^], yet we expect g[Fp^] < glFp'^], as 
it is easier to unpin from the shallower potential i. 

Now consider how the number of pinned vortices with thresholds in the range Fp < 
Fp^ < Fp + dFp changes during one time step. There are Ng{Fp)dFp such vortices at the 
start of the time step. According to rule two in §2.2[ Nfg{Fp)dFp vortices unpin from this 
threshold range (or indeed any other threshold range of width dFp, as rule two is unbiased 
with respect to Fp in this paper) due to thermal creep. According to rule three, all the 
remaining N{1 — f)g{Fp)dFp vortices unpin if Fp is less than Fm, an eventuality which occurs 
with probability dFuipiFu)', otherwise, none unpin. Finally, according to rule four, a 
number of vortices repin at sites in the threshold range [Fp, Fp + dFp]. The number is clearly 
proportional to A^0(Fp)(iFp, as 0(Fp) is the threshold distribution presented to an unpinned 
vortex by the nuclear lattice as it is about to repin. The constant of proportionality A is 
determined by normalization. (Recall that is constant in this paper, as we are interested 
in the glitch dynamics over decades.) In summary, we can write down the following master 
equation for g{Fp): 

^^d[N9{Fp)dFp] 



dt 

AN^{Fp)dFp - Nfg{Fp)dFp 



POO 

- N{1 - f)g{Fp)dFp / dFM HFm) • (3) 

JFp 

In the stationary state, the left-hand side of ([3]) vanishes, and we obtain 



giFp) = A0(Fp) 



/ + (!-/) / dFuHFu) 



-1 



(4) 



with A fixed by the normalization condition dFpg{Fp) = 1. An explicit analytic formula 
for g{Fp) is presented in equation flA3l) of Appendix |X] for the particular choices of 0(Fp) 
and iP{Fm) given by ([1]) and ([2]) respectively. With these choices, g{Fp) is zero outside the 
interval |Fp — Fo| < A, where no pinning sites are available. 
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2.4. Glitch sizes 

The time-averaged probabihty distribution function of glitch sizes s, denoted by h{s), 
follows directly from g{Fp). Let s = Au/u be the fractional rise in spin frequency during a 
glitch, as measured in pulsar timing experiments. For simplicity, we take s to be the number 
of vortices that unpin during a time step, divided by the total number of vortices (pinned 
and unpinned) in the star. [In reality, the contribution to Au/u from each vortex is also 
proportional to the distance it moves radial l y befo re repinning; see ^ rule three in §2.2[ 



the worked example in § J4.4l and lAlpar et al.l (119961 ).] If the global Magnus stress applied at 



that time step is Fm, then the resulting glitch size is given by 



s{Fu) = ef + e{l-f) I dF^giF^) . (5) 



The first term on the right-hand side of ([5]) was omitted by ISneppen fc NewmanI (119971 ) ; it 



is small (/ <^ 1). An analytic formula for s{Fu) is presented in equation (lA4p of Appendix 



lAl With (p{Fp) chosen according to ([T]), one obtains s = ef for all Fm < Fq — A and s = e 
for all Fm > Fo + A. 

Equation defines Fm implicitly (and uniquely) as a function of s in the interval 
|Fm — FqI < a. Hence the probability of getting a glitch with size in the range [s, s + ds] 
equals the probability that the global Magnus stress lies in the interval [Fm(s), Fu{s) + dFu], 
i.e. ip{FM)dFM_, which transforms into ip[Fu{s)][dFM{s) / ds]ds after changing variables by 
applying the chain rule. Combining this result with we arrive at 

= <(i - fHFu(s)] ■ 

Equation IQ can be evaluated formally by inverting ([5]) to obtain Fm(s) or, more practically, 
by evaluating s(Fm) from ([5]), h{Fu) from (Q, and then graphing h{s) parametrically. An 
analytic formula for h{s) is presented in equation (lASp of Appendix \M Care must be 
exercised wherever dFM{s)/ds diverges and s(Fm) is not invertible. For example, with 0(Fp) 
chosen according to ([1]), h{s) is bracketed by two delta-function spikes, (1 — e~^^''^^)5{s — ef) 
and e~^°~^5(s — e), corresponding to events with Fm < Fq — A (thermal creep only) and 
Fm > Fq -|- a (all vortices unpin) respectively. Even if 0(Fp), and hence g{Fp), are nonzero 
for all Fm [e.g. if 0(Fp) is Gaussian], the spikes remain: the delta functions are smeared out 
somewhat, but h{s) still diverges steeply (and integrably) as s ^ ef and s e. It is easy 
to verify this result analytically or with Monte-Carlo simulations. 
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Glitch statistics 



The glitch size distribution h{s) is a power law over many decades for / <C 1 (ISneppen fc Newman 



19971 ). Importantly, there are maximum and minimum cut-offs to the glitch size, which set 
natural limits on the extent of the power law. Its shape and normalization contain infor- 
mation about the pinning parameters in 0(-Fp). In this section, we demonstrate that it is 
possible in practice to infer these parameters by fitting the theory to observational data. In 
§3.11 and §3.2[ we implement the four rules in §2.21 in a simple Monte-Carlo simulation and 
study systematically how h{s) varies as a function of Fo/a and A/a. The simulation results 
are independent of for Nf ^1. In §3.31 and §3.4[ we fit the theoretical form of h{s) 
to data from the nine most active glitchers, both Poissonian and quasiperiodic, in order to 
derive constraints on Fo/a, A/cr, e, and /. 



3.1. Monte-Carlo simulations 

Figure [1] (left panel) displays a time series of the glitches generated by the automaton 
over the interval < At < 255 for = 10^ / = lO'^, e = lO'^, Fq = 4a, and A = O.GFq. 
The glitches occur intermittently, with a Poissonian waiting-time distribution, as arranged 
by construction through ([2]). Of the 255 events appearing in the left panel of Figured], most 
involve thermal unpinning only [Fu < F^ — A, s = ef = 10~^). However, there are 45 events 
where the Magnus force exceeds the pinning threshold at some of the occupied pinning sites 
(Fm > Fo — A), including one where the Magnus force unpins every vortex (Fm > Fq + A, 
s = e = 10-2). 

A frequency histogram of glitch sizes over a longer time interval (0 < At < 1 x 10®, 
corresponding to 10® events) is constructed for the same parameters and plotted on a log-log 
scale in Figure [T] (right panel, solid curve). The analytic form of h{s) derived from is 
plotted as a dashed curve over the simulation output. The two curves agree at the upper end 
of the s range, but are significantly discrepant at the lower end, because the time-averaged 
analytic theory in §2.31 does not capture the excess of small glitches arising from temporal 
correlations (aftershocks; see §1]). The discrepancies are smallest for A ~ Fq and Fq > a, 
the regime of interest when fitting to pulsar data. 

In the interval —4.5 < logio ^ ~ —2.5, h{s) is a power law with exponent —1.43 ± 0.01 
(weighted least squares fit). In addition, h{s) boasts two spikes at s = e/ and s = e, 
containing fractions 1 — e-(-^''-A)/cr ^ g gg g^^^^ ^~{Fo+A)/a ^ -^ j x 10"^ of the total number 

of glitches respectively, as discussed in §2.41 The power law, which does not encompass the 
spikes at the upper and lower end of the glitch size distribution, is generic: in the regime 
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Fp ^ — crln/, where forced unpinning dominates thermal creep, we have g{Fp) oc exp(Fp), 
s(Fm) oc expfFiu) for Fm ^ Fq — A, and hence h{s) oc exp[— 2Fm(s)] oc asymptotically 
( ISneppen fc Newmanll 19971 ). Physically, the power law emerges as a historical effect, because 
the system has memory. For example, a middling value of Fm (say, Fm = Fq) may trigger 
a middle-sized glitch (say, s ~ e/2), if g{Fp) is fairly flat (e.g. if every vortex unpinned 
during the previous iteration of the automaton). But the same middling value of Fm may 
trigger a tiny glitch instead (say, e/ < s ^ e/2), if g{Fp) is grossly depleted in the range 
Fq — A < Fp < Fq (e.g. following a chance sequence of iterations with Fm ~ Fq). The latter 
outcome is more probable than the former without preferring any particular glitch size, so 
h{s) scales as an inverse power of s over a large portion of its domain. 



3.2. Shape of h{s) 

The coherent noise model is completely specified by four parameters, which together fix 
the shape of h{s): Fq/ct, A/cr, e, and /. The minimum and maximum glitch sizes produced 
by the model are given by ef and e respectively. The roles of Fq/ct and A/cr are more subtle. 
In combination, the latter two parameters determine the extent of the scale-invariant portion 
of h{s), as well as its log-log slope. Figure [2] illustrates the various distributions that arise as 
we vary A/Fq in the range 0.2 < A/Fq < 1.0, given Fq = 0.25ct (left), Fq = 1.0a (middle), 
and Fq = 4. Oct (right). There are two panels for each value of Fq. The lower panel, which 
is drawn with a log-log scale, displays h{s) for A = 0.2Fo (dark grey), A = O.6F0 (medium 
grey), and A = I.OFq (light grey), showing the simulation output (solid histogram) and 
analytic result (dashed curve) from Appendix |X] superposed. The upper panel, which is 
drawn with a log-linear scale, displays iplFu) as a black histogram, together with g{Fp) for 
the minimum and maximum values of A/Fq (color coded as in the upper panel). We run 
each simulation to get 10^ events. The distributions in the upper panels are snapshots at 
the end of the run, binned in units of O.Olcr. The step in g{Fp) marks the last value of 
Fm sampled from '?/'(Fm) before the end of the run (followed by random repinning). The 
distributions in the lower panels are built up over the entire run, binned in units of 0.01 dex. 

Certain trends in the shape of h{s) are evident from Figure [2l (i) The distribution 
steepens as Fq/ct increases, while A/Fq is held constant, and vice versa, (ii) The probability 
density of the smallest glitches excluding the left-hand spike, viz. h{s e/~), increases as 
Fq/ct increases, while A/Fq is held constant, and vice versa. In contrast, the total probability 
in the left-hand spike increases as (Fq — A)/ct increases, (iii) To maximize the scale-invariant 
portion of h{s) and push the power-law exponent towards 2.0, we require Fq > ct and 
A ~ Fq. (iv) A gentle cusp appears in h{s) at small s for A > 0.6ct, rendering the cumulative 
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probability distribution J^j, ds' h[s') more convex. 

The above trends are related and easy to understand physically. To get a steep power 
law, the memory effect described earlier [i.e. gradual depletion of g{Fp) over part of its 
domain] must be free to take hold. Yet it cannot do so properly if the threshold distribution 
is narrow (which happens when A ^ Fq), or if large values of the Magnus force {Fu > -fo+^) 
occur frequently and reset g{Fp) before it can be depleted (which happens when Fq <^ a). 
The trends are also easy to understand analytically. For small s, where thermal creep 
dominates forced unpinning, equations (jlj)-(l6l) from the time-averaged theory imply that h{s) 
is flat, whereas the full simulation (which preserves temporal correlations) predicts h'{s) < 
in the limit s —>■ ef, reflecting the enhanced incidence of small aftershocks following a large 
event (see §1]). Elsewhere, in the regime Fp <^ — crln/, where forced unpinning dominates 
thermal creep, (!1])-(I6]) reduce to ( ]A7|) . as shown in Appendix |Al From ( ]A7p . it is clear that 
the turnover to a power law occurs for s > e(e^^ — 1)"^; as a corollary, a power-law portion 
of h{s) only develops if we have A > 0.35cr. This is confirmed by Figure [2l For example, 
you can see the turnover at log^^Q s ^ —2.6 (—4.1) in the dark (medium) grey curves in the 
right panel of Figure [2] and at log^g s ~ —2.8 in the light grey curve in the middle panel. The 
light grey curve in the right panel turns upwards at logj^g ^ ~ —2.8 because condition (1A6P 
on the smallness of /, which must be met to achieve a power law, is violated for / = 10"'^ 
and Fo = A = 4a. 



3.3. Extracting pinning parameters from observational data 

We are now equipped to fit the coherent noise model to glitch data from individual 
pulsars, in an effort to constrain the fundamental parameters F^/a, A/a, e, and /. By 
way of preparation, we note three things. First, A and hence a are directly measurable 
from waiting-time data. Second, the results in §3.11 imply that any observed size distri- 
bution is reproduced in shape by a relatively compact set of Fq and A values, which is 
encouraging. Third, the sizes of the smallest and largest glitches observed, denoted by 
(Az//z/)min and (Az//z/)max respectively, constrain e and / according to (Az//z/)max < e < 1 
and N-^ < f < (Az//z/)^in/(Az//z/) 

max- Of course, if 6 is much larger than (Az//z/)niax; 
we expect to see glitches larger than (Az//z/)max if we observe for longer. Luckily, this 
missing information is not a serious problem when fitting the model to the data, because 
I{Au/u) usually small. On the other hand, it is also possible that / is much 

smaller than (Az//z/)min/(Ai^/i^)inax, yet we fail to see glitches smaller than (Az//z/)min be- 
cause the resolution of pulsar timing experiments is imperfect. This problem is more serious, 
because J^^'^/'^^'^"^ ds' h{s') is usually large; much, perhaps most, of the distribution may be 
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invisible. To handle the proble m optimally, one should q uantify the resolution experimen- 



tally for each individual pulsar (iJanssen fc Stappersll2006l ). derive an observational window 



function w{s), and fit the data to h{s)w{s). Calibrated window functions are not published 
for most pulsars, so implementing the foregoing procedure lies outside the scope of this paper. 
Instead, we proceed conservatively by tailoring our fits to the observed bounds according to 
e = (Az//z/)max and / = (Ai//i/)min/(Az//z/)max- The results are interpreted critically in §4] 
in light of the uncertainties flagged above. 

In Figure [3], we fit the model to measurements of h{s) in the seve n pulsars which hav e 



glitched more than five times without a discernible periodic component (iMelatos et al.ll2008l ). 
i.e. with a purely Poissonian waiting-time distribution, as assumed by the model through 
([2]). The fits are done using a maximum likelihood approach, described in Appendix |Bl 
and the fitted parameters are recorded in Table [H There are two panels attached to each 
object in Figure [31 In the left panel, the likelihood function C is presented as a greyscale 
plot in the Fo/a-A/ Fq plane, such that larger C values are shaded darker than smaller C 
values. In the right panel, the continuous cumulative distribution J^^ ds' h{s') corresponding 
to the maximum likelihood — the best fit — is plotted as a dashed curve, together with the 
observational data (asterisks). It should be noted that the delta-function peaks that appear 
in the differential probability distributions at the upper and lower bounds on the glitch 
size appear as steps in the cumulative distributions. The objects in Figure |3] are ranked in 
decreasing order of the number of glitches observed, Ng. 

PSR J0534+2200 and PSR J1740-3015 are fitted well by the coherent noise model. 
For both objects, the best fit is achieved for Fq ^ 2a and A ^ O.QFq. The other objects 
are fitted reasonably well too, allowing for the simplicity of the model. For example, the 
fit for PSR J1825— 0935 (fifth row, left panel. Figure [3]) looks mediocre to the eye. But 
the likelihood function (right panel) is quite bumpy, with several peaks of almost equal 
heights, some of which (e.g. at A ^ 0.9-Fo) produce a fit which looks better to the eye, at 
the expense of a marginal drop in likelihood. For all seven objects, the best fits are achieved 
for A > 0.5Fo and 0.8 < Fq/ct < 3. This finding is in accord with the theoretical analysis in 
§3.11 one requires Fq > a and A ~ Fq in order to get a power-law size distribution h(s) oc s° 



that is steep enoug h (—2 ^ a < —1) to match the observations (IJanssen fc StappersI 12006 



Melatos et al.ll2008l ). In addition, one requires Fq + A > —a In / ^ {5—7)a in order to get a 
cumulative distribution f^j, ds' h{s') that is concave up. 

Some of the tightest constraints come from the spike in h{s) at s = ef. Four objects 
(PSR J0534+2200, PSR J1740-3015, PSR J1801-2304, and PSR J1825-0935) exhibit no 
evidence of a spike at s ^ (Az//z/)min in the data, implying either that these objects have 
A ^ Fq, or that (Ai^/i^)min comfortably exceeds ef because the experimental resolution 
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Fig. 1. — Sample output from a Monte-Carlo simulation of the coherent noise model with 
N = 10^ e = 10-2, / = 10-^ Fo = 4a, and A = O.GFq. Left panel. Time series of 
normalized glitch sizes, Au/u, over the time interval < At < 255. Right panel. Probability 
density function of normalized glitch sizes, p(Az//i^), accumulated over the time interval 
< At < 1 X 10^, showing the simulation output (solid histogram) and the analytic theory 
(dashed curve) overlaid. The spikes at Au/u = 10^^ (thermal creep only) and 10"^ (all 
vortices unpin) are real features, as explained in §2.41 
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Table 1: Pinning parameters extracted from the maximum likelihood fits in Figures [3] and 
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Fig. 2. — Theoretical size and pinning threshold distributions as functions of the pinning 
parameters Fq and A. Upper panels. Probability density functions of the Magnus stresses 
'^(Fm) (black histogram) and pinning thresholds g{Fp) (grey histogram), for A/Fq = 0.2 
(dark grey) and 1.0 (light grey), and for F^/a = 0.25 (left column), 1.0 (middle column), and 
4.0 (right column). Lower panels. Probability density function of glitch sizes for A/Fq = 0.2 
(black), 0.6 (dark grey), and 1.0 (light grey), and for Fq/ct = 0.25 (left column), 1.0 (middle 
column), and 4.0 (right column). The simulation output and analytic theory are graphed as 
solid histograms and dashed curves respectively. 
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Fig. 3. — Fits of the coherent noise model to pulsar ghtch data for the seven pulsars that have 
glitched at least six times, and whose waiting-time distributions do not show a discernible 
periodic component. Left column. Greyscale plots of the relative likelihood function for 
combinations of Fo/a and A/Fq in the ranges 0.1-5.0 and 0.5-1.0 respectively. The grey 
scale runs from least (white) to greatest (black) likelihood. The values of e and / used to 
generate the fits are printed on each panel, together with the maximum relative likelihood for 
the best fit. Right column. Cumulative probability function for the best fit (dashed curve), 
plotted over the measured cumulative probability function (asterisks). 
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prevents smaller glitches from being seen. On the other hand, the remaining three objects 
(PSR J0358+5413, PSR J0631+1036, and PSR J1341-6220) do exhibit some evidence of a 
spike at s ~ (Az//z/)min, implying A < O.6F0 and that the smallest glitches are resolved. By 
contrast, the spike at s = e does not constrain the fits tightly because it does not contain 
much integrated probability. 

The objects in Figure [3] and Table [1] are ranked in order of decreasing Ng to make the 
obvious but important point that the localization of the peak in the likelihood function, and 
hence the significance of the fit, improve markedly as A'^g increases. Encouragingly, the model 
also fits the data better as A^^g increases. The next task is to test the model more stringently 
by finding more glitches, both by reanalyzing timing data from the seven pulsars in Figure [3] 
to search for small, overlooked events, and by undertaking high-duty-cycle timing campaigns 
in the future. 



3.4. Quasiperiodic glitchers 



The timing histories of PSR J0537-6910 (A^g = 23) and PSR J0835-4510 (Vela; 
Ng = 17) harbor a periodi c glitching component, which coexists with the scale-invariant 



2006 



component discusse d so far (iLyne et al.l 119961 : [Marshall et al.l 120041 : iMiddleditch et al 
Melatos et al.ll2008l ). The periodic component comprises ~ 25 per cent of events and is 



characterized by narrowly peaked size and waiting-time distributions. When we attempt to 
fit the theoretical h{s) given by (1A5|) to the data from these two objects, following the same 
procedure as in Figure [3] and Appendix [Bl the agreement is no better or worse than for a 
pure Poissonian glitcher. Taken at face value, this is surprising; iIj{Fm), a basic input into 
the model, is wrongly specified by ([2]) for quasiperiodic glitchers, because it does not contain 
a narrowly peaked component. However, we show below that this omission does not show 
up noticeably in the observed h{s) for small A'g (e.g. A'g < 23). 

We recalculate the coherent noise model for V(-^m) = (1 ~ C'q)cr^^ exp(— Fm/ct) + 
Cqa~^S{FM/(J — Fq/a), where Fq is the M agnus force bu il t up during one period tq, and 



Cq defines the periodic fraction; see §5.2 of iMelatos et al.l (120081 ). Figure H] illustrates how 
the coherent noise process operates in the^ interesting case where F^ lies within the domain 



\Fp — Fq\ < A, where g{Fp) is nonzero. @ The model parameters are A^ 



10^ 



10 



-2 



For Fq < _Fo — A, the results in §3. II and §3.2l carrv over without change, because the periodic component 
does not unpin anything, and the number of ghtches with s > ef decreases by a factor Cq. For Fq > Fo + A, 
the periodic component unpins every vortex, and the number of ghtches with s = e increases by a factor Cq, 
while the rest of h{s) remains unchanged. 
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/ = 10-3, Fo = 4.0a, Fq = 4.0(t, and Cq = 0.25. The left panel of Figure 1 displays ?/'(Fm) 
on a log-linear scale; the spike housing the periodic component is clearly visible. The middle 
panel displays h{s) from the simulation on a log-log scale for A = 0.2Fq (dark grey his- 
togram), A = 0.6-Fo (medium grey histogram), and A = I.OFq (light grey histogram), with 
the analytic prediction for A = O.GFq overlaid as a dashed curve. The right panel displays 
J^j, ds' h{s') on a linear-log scale for the same three values of A (color coded as in the mid- 
dle panel), showing the simulation output including (solid curves) and excluding (dashed 
curves) the periodic component. One immediately sees that, for A = O.Qa and l.Ocr, where 
-Fq lies within the domain of g{Fp), there is a turnover in h{s) at s{Fq). Additionally, if 
the probability of getting Fm = Fq — A greatly exceeds that of getting Fm = Fq + A, such 
that thermal creep unpins more sites on average than events with Fm > Fq + A, then con- 
dition ( ]A6[) is violated, the power-law form of h{s) breaks down, and an excess (i.e. hump) 
of glitches emerges at small s. 

The main new feature in h{s) is a sequence of spikes at s = qe, g^e, q^e, . . . , with 
q = {Fq — Fq + A)/ (2 A), whose heights are in the ratio 1 : Cq : Cq : . . . . The spikes 
correspond to one, two, three, . . . consecutive draws from the periodic peak at Fm = Fq 
in ^(Fm), with each draw unpinning a fraction q of the vortices, of which a fraction q 
subsequently repin with Fp < Fq. The spikes in the simulation output do indeed occur at 
the predicted positions, although, apart from the first few, they are hard to see, swamped 
by the flood of Poissonian low-s events. On the other hand, the analytic, mean-field theory 
only predicts one spike, at s = qe; as a time-averaged theory (see §2.3p . it does not recognize 
correlated event sequences (e.g. consecutive Fm = Fq draws) as special. For example, the 
analytic theory does not distinguish between two Poisson events followed or separated by 
a periodic event, yet the relative frequency of these sequences determines the shape of the 
broad turnover in g{Fp) at Fp = Fq (where the analytic model predicts a simple step). 

The spikes are a key prediction of the coherent noise model for quasiperiodic glitchers. 
Do they show up in the data (as steps in the cumulative size distribution)? No; nor should 
we expect them to, when so few glitches have been detected, and each spike contains modest 
probability. To date, the objects PSR J0537— 6910 and PSR J0835— 4510 have been seen to 
sample the underlying event distribution A^^g = 23 and Ag = 17 times respectively. Monte- 
Carlo realizations of the coherent noise model with Cq = and Cq = 0.25 are statistically 
indistinguishable for A'g = 23 (and hence A'g = 17). Consequently, we experience no loss 
of generality if we fit a purely Poissonian model to the data from PSR J0537— 6910 and 
PSR J0835— 4510, following the same procedure as in FigureOand Appendix [Bl The results 
are graphed in Figure O (copying the format of Figure [3]), and the best fit parameters are 
recorded in the lower part of Table [TJ A respectable fit is achieved for PSR J0835— 4510, 
with Fq > a and A ^ Fq as usual, but not for PSR J0537— 6910, whose measured size 
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Fo/a=4,0 Fo/ct=4.0 Fo/ct=4.0 




Fig. 4. — Output of the coherent noise model including a periodic component weighted at 
25%, with N = 10^ e = lO'^, / = lO'^, Fq = 4.0a, Fq = 4.0a, and A/Fq = 0.2 (light grey), 
0.6 (dark grey), or 1.0 (black). Left panel. Magnus force probability density function, iplFu)- 
The spike at Fm = contains ~ 25% of the events. Middle panel. Glitch size probability 
density function h{s), featuring the output of the Monte-Carlo simulations (solid histogram) 
and the analytic theory for A = O.6F0 (dashed curve). Note the row of spikes generated 
by consecutive periodic events, whose s-positions and heights form geometric sequences 
with common ratios q and Cq respectively. Note also the hump at low s, which makes 
the cumulative distribution concave at low s. Right panel. Cumulative size distribution 
J^j,ds'h{s') for the models in the middle panel (matching color code). The solid (dashed) 
curves refer to models that include (exclude) the periodic component. 

distribution is skewed by two isolated glitches with s ~ 10~^'^ and a large group of glitches 
centered at s ~ 10~^'^ (implying, curiously, the existence of a larger periodic component than 
the observed waiting-time distribution can accommodate). The results raise the possibility 
that some small glitches were missed when reducing timing data from PSR J0537— 6910. 
As matters stand, the data are equally consistent with the presence or absence of periodic 
spikes. We need to detect more glitches with better resolution to settle the issue conclusively. 



4. Discussion 



In this paper, we propose an idealized model of pulsar glitches, i n which the collective 
unpin ning of superfluid vortices is treated as a coherent noise process ( ISneppen fc Newman 



19971 ). The model accounts for the scale invariance of glitch sizes over the range of possible 



glitch sizes, through the interplay between thermal and forced unpinning, with the two 
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Fig. 5. — Fits of the coherent noise model to ghtch data for PSR J0537— 6910 and PSR 
J0835— 4510, whose waiting-time distributions feature a periodic component. Lejt column. 
Grayscale plots of the relative likehhood function for combinations of F^jo and A/Fq in the 
ranges 0.1-5.0 and 0.5-1.0 respectively, in the context of a purely Poissonian model. The 
grey scale runs from least (white) to greatest (black) likehhood. The values of e and / used 
to generate the fits are printed on each panel, together with the maximum relative likelihood 
for the best fit. Right column. Cumulative probability function for the best fit (dashed 
curve), plotted over the measured cumulative probability function (asterisks). 



mechanisms balancing each other in a time-averaged sense. It reproduces approximately 
the measured size distributions of the nine most active glitchers for a range of physically 
sensible pinning threshold distributions, despite grossly oversimplifying the microphysics. 
More importantly, it demonstrates a fundamental and counterintuitive point of principle, 
which holds irrespective of the detailed microphysics: scale-invariant glitches can arise from 
homogeneous collective unpinning. Coherent noise is therefore a viable alternative to the 



-22- 



large family of inhomogeneous mechanisms proposed in the hterature (e.g. self-organized 
criticality) , which rely upon nearest-neighbor avalanches and the existence of large-scale 
capacitive regions. 

Such an alternative is welcome. It fits with comtemporary notions of nuclear pinning and 
thermal creep, where it is thought that vortices hop between adjacent (or nearly adjacent) 
lattice pinnning sites, which lie ~ 30 fm apart and are therefore distribut ed homogeneously 



on relevant macroscopic sca l es, e.g. the perio d of the Abrikosov lattice (lAlpar et al.lll984 



Joneslll99ll : iLink et al.lll993l : |Jahan-Mirill2006l ). It also helps to explain the puzzle, posed in 
^ of why glitch sizes vary so much in an individual pulsar, even though the global Magnus 
stress builds up to a similar level between glitches (dictated by the Poissonian waiting-time 
distribution) and is felt simultaneously by every pinned vortex. 



4.1. Fn and A 



Disaggregated glitch statistics for individual pulsars became available in r neaningful vol- 
ume only recently, yet already they tell a clearer story than aggregate statistics (IMelatos et al. 
20081 ). (i) When disaggregated, the waiting-time distribution is Poissonian, except in two 
pulsars which have a small (~ 25%) periodic component, (ii) When disaggregated, the size 
distribution is consistent with a power law h{s) oc s", but the power-law exponent a is not 
universal; the aggregate size distribution is inconsiste nt with a surn of id entically sloped 
power laws at the 99% level of confidence [see §4.4 of iMelatos et al.l (120081 )]. (iii) The ob- 
served cumulative size distribution of the Poissonian glitchers is concave in two objects, 
convex in two others, and contains a point of inflexion in the remaining three. 

The glitch model in this paper is motivated by (i). Its predictions are consistent with 
(ii) and (iii). We find the following properties. First, the theoretical h{s) is scale invariant 
over several decades in each object, with —2 < a < 0. Second, the observed absence of a 
step at s = e/ in most objects implies a broad pinning threshold distribution, with A Fq. 
Physically, this conclusion is natural: nuclear structure calculations in the local density and 
Hartree-Fock-Bogoliubov approximations independently indicate that the pinning energies 
in a neutron star covers a br oad range (1 MeV < E-^ < 4Mey), because the superfiuid 
densi t y changes with depth (ILink fc Epsteinlll99ll : |Joneslll997l . Il998l : iDonati fc Pizzochero 



2003 



20041 . l2006l : lAvogadro et al.l 120071 ). translating into a broad 4>{F^) when we fit our 

homogeneous model to data. [f| Third, most objects have —2 ^ a ^ —1, i.e. h[s) is not 



^ When vortices repin after a glitch, g{Fp) is broad. During thermal creep, which involves continuous 
angular momentum transfer from the superfiuid to the crust, the vortices gradually bend in a piecewise 
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flat, which requires Fq > a. Fourth, objects with a concave segment in h{s) fall in the 
regime where the rates of thermal and forced unpinning are comparable [i.e. condition flA6p 
is violated], producing a hump in h{s) at small s. 

The coherent noise model also works fairly well for the quasiperiodic glitchers PSR 
J0537— 6910 and PSR J0835— 4510 and leads to similar conclusions, i.e. Fq > a and A ^ Fq. 
In principle, if ^(-Fm) contains a periodic component, a sequence of equally spaced (in log^^g ^) 
spikes appears in h{s) (see §3.41) . The available data are consistent with the spikes, as Figure 
O shows, but one cannot hope to actually trace out the associated steps in the cumulative 
size distribution unless A^g > 10^. The periodic component impacts more on the waiting- 
time distribution, which flattens at small At, than on h{s), which rises at small s, swamping 
the spikes. Crucially, quasiperiodic glitching does not discriminate between homogenoeus 
(e.g. coherent noise) and inhomogeneous (e.g. avalanche) collective unpinning; it is equally 
at home in both families of models. For example, in self-organized critical systems (e.g. sand 
piles), quasiperiodic i ty is caused by system-spanning av alanches triggered by a fast external 
driver JjensenlllQQsl : ISornettel[20o3 : iMelatos et allboosh . 



In order to exploit the results presented in Tabled] to infer an actual value of Fq, we must 
define a in terms of physical pulsar parameters. We present two possible definitions of a. 
Firstly, we assume that the Magnus stress that accumulates between glitches is proportional 
to the time between glitches, such that 



(Ti = 27Ti'RpsK,/ X 



(7) 



where A is the mean glitching rate, z> is the rate of angular deceleration, R is the char- 
acteristic radius at which vortex pinning occurs (or, more correctly, the characteristic ra- 
dius w here pinning occurs rnost s t rongly), and Ps is the superfluid mass density at that 
radius (IDonati fc Pizzocherd l2006l : lAvogadro et al.l 120071 ) . The Fm — At correlation un- 
derlying Equation ([7]) h as been proposed often in the literature in the context of models 
without thermal creep (lAnderson fc Itohl 1 19751 : iMorley fc Schmidtl Il996l : iLyne et al.l l2000l : 
Wang et al.l 120001) and f inds support in the reservoir effect observed in PSR J0537— 6910 
( iMiddleditch et al.ll2006l ). This is an opportune time to test its consequences. 



A key corollary of [7] is that the average pinning strength Fq ~ a varies greatly between 
neutron stars. In the left panel of Figure O we plot Fq (in Nm~^) versus characteristic 
spin-down age = — z//(2z>) (in kyr) for the nine pulsars in Table [T] assuming a = ai. We 



fashion, via a sequence of small rearrangements which bring them into line with local microcrystalline 
boundaries. Consequently, {Fp) incre ases a n d qiF „) narrows, until the next glitch intervenes. For a fuller 
exposition of this idea, please consult I Joned (jl99lh . 



-24- 



estimate the error bars conservatively to be ±0.5 dex from Figures [3] and [5] by inspection. 
Remarkably, there is an inverse correlation over three decades in and two decades in 
Fq, with Fq oc 7--i-28±o.3i approximately. Taken at face value, the correlation is surprising 
physically: Fq is determined by the structure of the nuclear lattice, averaged over position 
within the stellar crust in our model. It is therefore primarily a function of the mean crustal 
density, which varies slightly across the neutron star population. One might argue that 
Figure [6] indicates a temperature effect: as the star cools, thermal creep occurs more slowly. 
However, the rate of thermal creep is parametrized by /, not Fq. One might argue instead 
that Fo depends sensitively on and hence T. However, this effect seems to work in the 
wrong direction; iDonati fc Pizzochero found that Ep increases as T decreases and ps 

increases. We therefore speculate that the nuclear lattice "anneals" as the star ages, with 
weak seismic and/or thermal fluctuations erasing defects as time passes. This process of 
annealing would have a less dramatic effect on the mean pinning energy if pinning is mainly 
due to sites in the nuclear lattice. It is well known that the concentration of monovacancies 



affects the strength of pinning severely (|Joneslll99ll . 119971 . Il998l ). But it is an open question 
if a realistic annealing process exists to lower the monovacancy concentration with time. 
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Fig. 6. — Best fit mean pinning threshold Fq (units: Nm~^) versus characteristic spin-down 
age Tc (units: kyr) for all nine pulsars that have glitched at least six times for two different 
definitions of cr: di from Equation (left) and cr2 from Equation ([H]) (right). The data are 
drawn from Table [TJ Quasiperiodic glitchers are plotted in grey. 



An alternative definition of a invokes the process of ther mal creep to define the mean dif- 



feren t ial lag between th e interior superfluid and stellar crust f Alpar et al.lll984j ; lLink Sz Epstein 
199ll : Ijahan-Mirilbooeh . viz. 

(J2 = 6ujRpsH, (8) 
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with 

5uj = LUcr—T- In (AtcVo/R) , (9) 

Ep 

where tOcr ~ O.lrads^^ is the critical lag above which vortex pinning can no longer be 
sustained, ksT/Ep ^ 30 is the ratio of the typical thermal energy and pinning energies, 
and Wo ~ 10^ ms^^ is a trial microscopic vortex creep speed. [Note that no self-consistent 
theory exists at present for the stochastic fluctuations of 6uj from glitch to glitch about its 
mean Q.] When the fitted values of Fq are recalculated using a = cr2 and plotted against 
Tc (right panel of Figure [6]), the inverse correlation observed in the left panel of Figure [6] 
is no longer present at a significant level. In this case Fq oc r~^'^^^^'^^ , suggesting that the 
mean vortex pinning strength does not change from pulsar to pulsar and therefore seems to 
be independent of pulsar temperature. It should be noted, however, that both cxi and a2 are 
interpreted assuming that the characteristic radius at which pinning occurs does not depend 
on characteristic age (and hence temperature) and therefore it does not vary from pulsar to 
pulsar. 



4.2. e and / 



Figure [7] plots the lower and upper bounds implied by the data on e (left panel) and / 
(right panel) respectively against spin-down age Tc for the nine pulsars in Table [TJ There is 
no obvious trend, although a real trend can easily be concealed, if the bounds lie far from 
the true physical values for the resolution-related reasons outlined in §3.31 In particular, one 
expects / to decrease steeply as Tc increases. As a pulsar ages, it cools, and the thermal 
unpinning rate should fall exponentially. On the other hand, Fq also decreases with (Figure 
E]), compensating for cooling at least partly. 

We can estimate / crudely from the general Arrhenius formula for thermally activated 
processes and compare with the upper bounds in the right panel of Figure [71 Let us take 
/ ~ exp[ 
ujrr is the 



( Jones 


1991; 


Alpar et al. 


1996; 


Jahan-Miri 


2006) 



Ep^2 MeV and T 



10^ K, we obtain Ep/k-QT ^ 2x 10^. This implies a creep rate / which 
is certainly consistent with, but far below, the upper bounds / < 10~^ in the right panel of 
FigureO To approach / ~ 10~^, one requires either |Fo— A| ~ 10~'^Fo (/ is dominated by the 
shallowest pinning sites) or l — u/ucr ~ 10""^ (most of the stored differential rotation persists 
after a glitch). The upper bounds approximate the true value closely, if the experimental 
resolution is much better than (Az// z/)min [and e ~ (Az// z/)max], because smaller ghtches occur 
more commonly than larger glitches and are therefore certain to be seen if the ex perimental 



resolution allows. For example, in PSR J1740— 3015, iJanssen fc Stapperj (120061 ) simulated 
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Fig. 7. — Left panel. Minimum fraction of pinned vortices as a proportion of the total, 
e < (Az//i/)max, plotted versus spin-down age, measured in kyr. Right panel. Maximum 
fraction of pinned vortices that unpin thermally every time step in the coherent noise model, 
/ > (Az//z/)„iin/(A^/'^)max) plottcd vcrsus spiu-dowu age, measured in kyr. 



microglitch detection in a noisy signal and estimated the resolution for that object to be 
1 X 10~^^, well below {Au/u)ram = 7 x 10~^°. Microglitch detection simulations are needed 
for other objects. § 

It should be noted that vortex creep and inhomogeneous unpinning are not independent 
processes. In fact, a continuous-time coherent noise model, in which vortices unpin thermally 
at a continuous rate given by the Arrhenius formula, has the potential to include both 
phenomena self-consistently as opposite extremes of the unpinning dynamics. That the 
current model does not model thermal unpinning explicitly is a direct result of it being 
discrete in time. 

The data easily accommodate e ^ (Ai//i^)max, because ds' h{s') is typically small for 
s > (Ai//t/)nia,y H owever, it is debatable whether this is actually necessary. On the one hand, 
Lyne et al. measured e = 0.017 ± 0. 002, well abov e the l argest glitch ever observed 



in any pulsar [(Ai//i/)max = 2 x 10~^; see iMelatos et al.l (120081 )]. One can also argue for 
e ~ 1 on physical grounds in the context of the coherent noise model, where lattice sites and 
defects are microscopically separated and m uch more numerou s than vortices. On the other 
hand, the aggregate value of e measured by iLyne et al.l (120001 ) effectively averages together 
different pulsars, binned over semi-decades in z>. While this approach reduces the formal 
error bar, it obscures the physical interpretation, given the likelihood that e differs from 



Alpar fc Bavkall ( 20061 ) argued that anomalous brakin g: indices observed in several pulsars are partly 



attributable to unobserved microglitches. Jahan-Miril ( 20061 ) related the creep rate to the time- averaged u. 
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pulsar to pulsar, and that (Az/) is dominated by the largest, rarest, and hence unobserve d 
(over 40 yr of monitoring) glitches; see the detailed discussion in §6 of iMelatos et al.l (120081 ) . 
It is currently possible to measure e reliably in only one object, PSR J0358+5413, which has 
e < 7 X 10-5 jMelatos et al.lboosh . 



4.3. Aftershocks: an observational test of homogeneity 

To the eye, the glitch statistics in Figures [1] and [2] are indist inguishable from those 



produ ced by a nearest-neighbor avalanche process, e.g. Figure 1 in IWarszawski fc Melatos 



( 120081 ). Of course, the underlying collective behavior is very different: in coherent noise, the 
spatial correlation function Gc{i,j) = (Fp^Fp^) — (Fp^)'^ is independent of vortex separation 
|xj — Xjl, whereas a self-organized critical system with nearest-neighbor avalanches exhibits 
long-range correlations on all scales, with Gc oc |xj — Xjl"'^ and f3 > 0. Unfortunately, 
we cannot measure Gc directly in a neutron star. However, it turns out that we can still 
discriminate between homogeneous (e.g. coherent noise) and inhomogeneous (e.g. avalanches) 
glitch mechanisms observationally, by taking advantage of temporal correlations embedded 
in the data. Below, we propose three related observational tests of this kind, all of which 
are practical in the medium term. 

1. Aftershocks. Aftershocks occur in a coherent noise process because, after a large glitch 
with Fm > -Fo + A, the threshold distribution g{Fp) is repopulated evenly; that is, {Fp) 
immediately after a large glitch is less than (Fp) for the time- averaged g{Fp), which 
is depleted at l ow Fp. Consequently, more vortices than usual unpin during the next 



few time steps (ISneppen fc NewmanI 119971 ). In contrast, aftershocks are not produced 



by nearest-neighbor avalanches, because successive glitches arise from the relaxation 
of insular capacitive domains and are therefore independent. 

Aftershocks cannot be analyzed by the time-averaged, mean-field, analytic theory in 
§2.3^ which contains no information about temporal correlations. U But Monte-Carlo 



"^As defined in lAlpar et all 1 19961 . and references therein), capacitive domains may indeed be connected, 



and hence their relaxation not independent. In the wider (e.g. experimental) SOC literature, capaci- 
tive domains are def in ed such that all connected "subdomains" constitute one domain. See for example 
Sornette et al. ( 1991 ): JensenI |"l998 l. Defined this way, the relaxation of successive domains is statisti- 



cally independent (exce pt for system- wide avalanches), as verified by experiments with sand and rice grains 
iRosendahl et al.lll993h . 



^ The time-averaged theory assumes a unique glitch size s{Fm) for every Magnus force Fm, whereas 
in reality the automaton in i )2.2l produces a power-law distribution of sizes for fixed -Fm; see Figure 8 of 
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simulations of the automaton in §2.21 reveal the aftershock effect clearly. Figure [8] 
displays the conditional size probability density function h{s,t + At\s',t) for s' > e, 
0.67e, 0.5e, and 0.33e (color coded from dark to light) and compares it with the full 
size distribution h{s) = h{s,t + At\s' ,t) (lightest grey). Clearly, there is a relative 
excess of large glitches following a large glitch, especially for s' = e. Note that the 
aftershocks discussed here are standard glitches; they are not the sarn e as th e time- 
resolved secondary spin up events noted in the Crab by IWong et al.l (120011 ). which 
occur 20-40 d after a glitch. 
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Fig. 8. — Conditional size probability density function h{s, t + At\s', t) for s' > e (thick black 
histogram), 0.67e (dark grey), 0.5e (medium grey), and 0.33e (light grey), together with the 
full probability density function h{s) = J2s' h{s,t + At\s' ,t) (lightest grey). The conditional 
function describes the probability of getting a glitch of size s at time t + At following a glitch 
of size s' at time t. The excess of large glitches following a large glitch is a clear signature 
of the aftershock effect. 



2. Glitch lifetimes. The distribution of glitch lifetimes, or rise times (i.e. durations; cf. 
waiting times), is also a power law i n a coherent noise process , with exponent ^1.0 
(plus a small logarithmic correction) (jSneppen fc Newmanlll997l ). The power law stems 



Sneppen fc Ne'wman 



(|l997l ). 
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from the memory effect described in §3.11 Its exponent is the same for all forms of 
ijj{Fy[)- In c ontrast, an avalanche process leads to a range of exponents for the lifetime 
power law (IWarszawski fc Melatod l2008l ). Testing this prediction is ambitious. At 
present, the best measurement of glitch d urations is ari upper limit of ~ 40 s, which 
comes from continuous monitoring of Vela (jDodsorul2002l ). However, it may be possible 
to improve on the resolution currently available through single pulse timing and daily 
monitoring with low-radio- frequency arrays currently un der construction [M. Bailes 
and R. Bhat, private communication, iDodson et al.l (j2009l )]. 

It should be noted, however, that a single component model like the one presented here 
cannot accurately describe glitch lifetimes, which depend on the coupling between 
the stellar interior and its crust. This coupling is thought to rely on the presence 
of superfluid neutrons and superconducting protons in the core of the neutron star 
( lAlpar fc Baykalll2006l ). A two component coherent noise model, capable of addressing 
the issue of crust-core coupling, will be the subject of a future paper. 

3. Size versus waiting time. The coherent noise model predicts a statistical correlation 
between the size of a glitch and the time since the previous glitch. In contrast, an 
avalanche process exhibits no such correlation. For example, it can store Magnus stress 
in metastable reservoirs over a long period of time, punct uated by min o r glitches , until 



a major glitch releases most of the stress all at once (jJensenl Il998l : ISornettd 12004 



Melatos et al.l l2008l : IWarszawski fc Melatod 120081 ) . The necessary homogeneity of the 



model presented here does not allow for the unique treatment of different c omponents 
of the stellar interior. If, as is discussed in detail in lAlpar fc Baykall (120061 ) . we could 
differentiate between the capacitive and resistive regions of the superfluid, we would 
expect similar correlations between the pulsar spin-down rate, the interglitch waiting 
time and glitch size. In particular, by identifying regions of the superfluid where 
fluctuations in vortex pinning arise thermally and regions where vor tices are always 
pinne d, the effective pinned fraction e becomes spatially dependent (lAlpar fc Baykal 



20061 ). 



The coherent noise correlation is imperfect because of the memory effect described in 
§3.1j two events with the same Fm can lead to different s values, according to whether or 
not g{Fp) is depleted at low s by a prior sequence of relatively small Fm. Nevertheless, 
Figure M demonstrates that it is present on average. In the upper panel of Figure [HI 
the size of a glitch is plotted against the time since the previous glitch (normalized 
to A^^) for a representative Monte-Carlo simulation of the coherent noise model. The 
associated linear Pearson correlation coefficient is r = 0.57 for 10^ events, implying a 
substantial correlation. In the lower panel of Figure [HI the same plot is constructed 
for the 23 glitches of PSR J0534+2200 (black asterisks). The associated Pearson 
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coefficient is r = —0.22, implying a weaker correlation. Likewise, we find r = 0.16 for 
PSR J0835— 4510, which is not plotted. We advise caution when interpreting these 
results, (i) For small A^g, a sample of A^g events does not exhibit the correlation in 
the upper panel at a statistically significant level, (ii) The correlation stems from 
the aftershock effect. It is dominated by the small-s events populating the dark band 
in the upper panel, the part of the distribution that is measured least reliably (e.g. 
unobserved microglitches). 



Wang et al. 


(2000 


) and 


Wone et al. 


(2001) 



However, all previous studies have analyzed the statistics from all glitching pulsars in 
aggregate (grey asterisks in the lower panel of Figure [9]) , in an e ffort to maxiniize th e 



number of data points. Since pulsars have different A and (Az/) ( iMelatos et al.ll2008l ). 
both axes of the potential correlation are washed out in aggregate data. It is therefore 
vital to test pulsars individually. In the case of the much-studied Vela pulsar, models 
that include correlations between the spin-down rate and the interg litch waiting time 
have successfully describ ed ghtch behav i our (e .g. lAlpar et al.lll984). S i milar analyses 



have been conducted by [Marshall et al.l (12004 ) and iMiddleditch et all (l2006f ) on PSR 
J0537-6910, resulting in claims that the interglitch interval can be predicted to within 
a few days. 



For the definition of a that assumes that the Magnus stress is proportional to the time 
elapsed (o"i), the coherent noise model analyzed in this paper is incomplete in one major 
respect. It does not predict endogenously the waiting-time distribution and hence the mean 
rate A. Rather, A is measured observationally and put into the model by hand through 
ipi^Fyi). There is nothing wrong with this approach, of course; it takes advantage of a well- 
determined observational fact to construct the theory. But it does mean that the model is 
incomplete. In particular, as it stands, the model cannot answer important questions like 
why some pulsars glitch and others do not, and why, of the objects that do glitch, some are 
more active than others. Moreover, it can only be applied to objects where enough glitches 
have been observed to measure A reliably. On the positive side of the ledger, there is no 
reason in principle why a microphysical theory of the waiting-time distribution and hence A 
cannot be developed by studying the temporal behavior of single- vortex unpinning triggers 



like thermal activation (lAlpar et al.l Il984 



Joned Il99ll : iJahan-Miril |2006| ) . Once available 



such a theory, combined with the coherent noise process to explain collective unpinning, and 
generalized to include the radial dependence of 0(-Fp), may offer a path to a complete glitch 
theory. The beginnings of such an approach are contained in Equations ([8]) and ([9]) in § 14.11 
A similar approach may also prove fr uitful in interpret i ng classic experiments on magnetic 
flux creep in type II superconductors (iField et al.lll995l : iBassler &: Paczuskilll998l ). 
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Fig. 9. — Upper panel. Glitch size, s, versus the normahzed time since the previous ghtch, 
AAt, for a representative simulation of the coherent noise model, with N = 10^, e = 10~^, 
/ = 10~^, Fo = 4.0(7, and A = O.GFq. Each dot marks one glitch. The distribution is 
bounded above by the function ef + e(l — f)[XaAt — {Fq — A)]/ (2A) (dashed curve) and is 
centered around the time-averaged relation s(Fm = X/aAt) from equation ([5]) (dotted curve). 
The Pearson linear correlation coefficient r is printed on the plot. Lower panel. Glitch size, 
s, versus the normalized time since the previous glitch, XAt, for PSR J0534+2200 (black 
asterisks) and for every glitch observed to date (grey asterisks). 

4.4. Observational verification 

The coherent noise model in this paper is unrealistic in three important respects, which 
render it difficult to verify observationally. (i) The model is discrete, not continuous, in time, 
(ii) It is homogeneous, so it cannot include the radial dependence of pinning strength (or creep 
rate), and the moments of inertia of the regions from/through/to which the vortices move, 
in a self-consistent way. For example, vortex motion per se is not modelled; we do not track 
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Fig. 10. — Theoretical size and pinning threshold distributions as functions of the pinning 
parameters Fq and A. Parameters are the same as for Figure [21 but it is now assumed 
that unpinned vortices traverse at random one of four regions moments of inertia accounting 
for 0.15, 0.2, 0.25 and 0.4 of the total before repinning. The fractional glitch size is also 
multiplied by the ratio of the total stellar moment of inertia to the moment of inertia of the 
crust (/o//c)- 



where individual vortices go when they unpin, (iii) It allows no feedback of the vortex motion 
on the two stellar components (pinned superfiuid and crust), and hence does not determine 
the observable rotational dynamics in an internally self-consistent way. These weaknesses 
do not interfere with the main physical result: namely, and counterintuitively, that you get 
a power-law distribution of event sizes from a Poisson distribution of waiting times, in spite 
of the system being homogeneous, and without nearest-neighbor avalanches occurring. These 
three weaknesses significantly reduce the ease with which this mod el can be distinguished 



from other glitch models, such as those invoking avalanche dynamics (IWarszawski fc Melatos 



20081 . for example). The appearance of features such as aftershocks, correlations between 
glitch sizes and waiting times, and particular power-law exponents, may depend strongly on 
the coupling of the different stellar components, and hence require that the model account 
for this coupling. 

In a first, crude attempt to generalize the coherent noise model for multiple components, 
we augment the model automaton described in § 12.21 so that the vortices unpinned during 
each glitch to pass randomly through one of four regions of the star before repinning. The 
choice of fou r regions is a r bitrar y and merely illustrative. Following the vortex creep model 



presented in lAlpar et al.l (119861 ) . we assume that the size of a glitch depends on both the 



number of vortices that unpin, and the fraction of the total stellar moment of inertia through 
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which the unpinned vortices move before repinning. A homogeneous model does not record 
the starting and end points of the vortex motion during a ghtch, and hence cannot accurately 
predict the change in angular momentum of the superfluid resulting from the (un)pinning 
process. As a worked example, we assume that the fractional volume traversed by unpinned 
vortices falls into one of four broad bins, accounting for 0.15, 0.2, 0.25, and 0.4 of the total 
moment of inertia of the star Jq. The glitch size is then calculated as the product of the 
unpinned vortex fraction and traversed moment-of-inertia fraction. The results of Monte- 
Carlo simulations using the modified "multi-component" automaton are shown in Figure [TUl 
The main effect, when compared to the single component model shown in Figure [2] for the 
same values of -Fq/ct, A/a, e, and /, is to blur the upper and lower edges of the power law 
(without changing much the exponent in between). That is to say, multiple components 
change the effective value of e. 

We note that in the realistic case where all unpinned vortices traverse similar volumes, 
the multiple spikes appearing in the size distributions in Figure [TOl merge into single peaks at 
both ends of the distribution. A self-consistent study of this important effect awaits future 
work. 

In closing, we emphasize again that the model in this paper is highly idealized. Quan- 
titative conclusions drawn from fitting the model to observational data should be viewed 
cautiously until more glitches have been observed. 

The authors thank Stuart Wyithe for expert advice on how to rigorously fit the theory 
to the data in ^ We thank the referee, Dr Ali Alpar, for his thorough reading of the paper; 
his insightful comments have greatly improved this manuscript. This research was supported 
by an Australian Postgraduate Award, the University of Melbourne Research Grant Scheme, 
and the University of Melbourne-CSIRO Collaborative Research Support Scheme. 
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A. Analytic results for a top-hat distribution of pinning thresholds 

In this appendix, we write down explicit analytic formulas for the main results of the 
mean-field theory when applied to a top-hat distribution of pinning thresholds, given by ([I]), 
and an exponential distribution of Magnus stresses, given by ([2]). To assist readability, we 
redefine the force (per unit length) variables Fm, Fp, Fq, and A, appearing in flA2l) - flA7l) 
below, to be normalized versions of their counterparts in the main text, measured in units 
of a. We also introduce the auxiliary function 



X{x) = 1 



/ + /exp(x) 



(Al) 



and the parameter 



/i = /i(/,Fo,A) 



In 



A(Fo + A) 
A(Fo - A) 



(A2) 



1. Time-averaged threshold distribution at occupied pinning sites: 




/exp(Fp)[^A(Fp)]-i 

xH{Fp -Fo + A)H{-Fp + Fq + A) . 



(A3) 



2. Glitch size as a function of Magnus stress: 



s{Fm) 




xH{Fm -Fo + A)H{~Fm + Fo + A) 
+efH{-FM + Fo - A) 
+eH{Fu - Fo - A) . 



(A4) 



This preprint was prepared with the AAS IM^t;X macros v5.2. 
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3. Time-averaged probability distribution function of glitch sizes: 



h{s) 



ffxXjFo - A) 



exp 



- ef) 
e(l - /) 



x<-(l-/) + A(Fo-A)exp 

xH{s-ef)H{-s + e) 
+ [l-exp(-Fo + A)]5(s-e/) 
+ exp(-Fo- A)5(s-e) . 



- ef) 
e(l - /) 



-2 



(A5) 



If / is small, such that the interval f < s/e < 1 spans several decades, h{s) tends to a 
universal power law, namely h{s) oc s~^, in the regime s ^ ef. To see this mathematically, 
consider the limit where / is small enough such that fi is also small. A sufficient condition 
for this ordering to obtain is 

/exp(Fo + A)<l. (A6) 
When (IA6I) is satisfied, we find fi ~ /[exp(Fo + A) — exp(Fo — A)] and hence 



his) 



eexp(-Fo + A) 
exp (2 A) - 1 



s + 



exp(2A) - 1_ 



-2 



(A7) 



for ef <^ s < e. A similar result follows for any other physically reasonable choices of (f){Fp) 
and V'(-Fm) ( jSneppen fc NewmanI 119971 ) . 



B. Maximum likelihood algorithm for fitting h{s) 

In this appendix, we describe briefly the fltting algorithm used to produce Figures [3] and 
|5] and Table [H to help the reader reproduce the results. The algorithm can handle relatively 
small data sets meaningfully, but of course the statistical significance of its output improves 
as Ng increases. 

Figure [11] demonstrates how to obtain a best fit in three steps, taking as an example the 
Crab (PSR J0534+2200), which has glitched A'^g = 23 times, (i) Starting with a particular 
choice of F^/a and A/Fq, we create many (~ 10^) realizations of the model by sampling the 
continuous, theoretical h{s) Ng times to construct each realization. A representative realiza- 
tion is plotted as a cumulative probability distribution in the left panel of Figure [H] (black 
asterisks), together with the continuous distribution j^^ds' h[s') from which the realization 
is sampled (dashed curve), and the observational data (grey asterisks), (ii) We compute the 
maximum unsigned separation D between each realization and the continuous distribution. 
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i.e. the Kolmogorov-Smirnov statistic, □ and bin the results to produce a frequency histogram 
of D which is unweighted by bin width and normahzed to have unit area, as in the middle 
panel of Figure [TTl (Note that the plotted histogram is representative; it does not corre- 
spond to the choice of Fq and A that gives the best fit.) We then compute the maximum 
unsigned separation -Ddata between the observational data and the continuous distribution. 
The height of the histogram oX D = -Ddata, where the dashed lines intersect in the middle 
panel of Figure [TT| is called the relative likelihood £(Fo/cr, A/Fq) for that particular choice 
of -fo/o" and A/Fq. (iii) We repeat for a range of F^/a and A/Fq. The best fit parameters 
are those that maximize C{FQ/a, A/Fq). For example, in Figure [TT| the maximum C = 0.23 
is achieved for {Fo/a, A/ Fq) = (1.8,0.92) (here, the gridding differs slightly from Figure [3]). 
The corresponding continuous cumulative distribution is plotted as a dashed curve in the 
right panel of Figure [TTl together with the observational data (grey asterisks). 



^ By tailoring ef to match {Ai'/i')^in, we artificially impose the restriction D > N^^, caused by the 
guaranteed discrepancy in the leftmost bin. Although the restriction varies from pulsar to pulsar, it applies 
equally to all the model realizations and the observational data in any individual object, so its distorting 
influence is mild. 

The relative likelihood C only has meaning when used to compare models with different (Fq/ct, A/Fq) 
in the same pulsar. It cannot be used to compare models across different pulsars. C is neither a probability 
nor a probability density. If we have C{Fo,D) = 1.5C{Fq, D'), say, we can conclude that model {Fo,D) is 
more likely than model {Fq,D'), but not that it is 1.5 times more likely. 
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Fig. 11. — Graphical demonstration of the maximum hkehhood algorithm for fitting 
Lejt 'panel. A model is chosen by selecting Fq and A. Many (~ 10^) realizations of the model 
are created by sampling the theoretical distribution iVg times. One representative realiza- 
tion is plotted here (black asterisks), together with the continuous cumulative probability 
distribution from which it is drawn (dashed curve) and the observational data (grey aster- 
isks). Middle panel. The Kolmogorov-Smirnov D statistic quantifies the separation between 
the realization and the underlying continuous distribution. A frequency histogram of the D 
statistic is constructed from all the reahzations (solid staircase), normalized to unit area. 
The separation between the data and the continuous distribution, Ddata, is also computed. 
The relative likelihood of a trial pair (Fq, A), denoted by C, is defined as the height of the 
histogram at D = -Ddataj where the dashed horizontal and vertical lines meet. Right panel. 
The foregoing procedure is repeated for many combinations of Fq and A, each time yielding 
a D histogram, a value of D^ata, and a relative likelihood C The best fit parameters max- 
imize jC. Here, the best fit is achieved for e = 2.1 x 10"''', / = 3.2 x 10~^, Fq/ct = 1.8, and 
A/Fq — 0.92. The associated cumulative probabihty distribution (dashed curve) is plotted 
over the observational data (grey asterisks). 



